Multi-objective optimization method for power system peak-shaving based on index linkage analysis

ABSTRACT

The invention relates to a multi-objective optimization method for power system peak-shaving based on index linkage analysis. The invention comprehensively considers five peak shaving indexes, and introduces a linkage analysis theory in economics to test stationarity and causality of data series of peak-shaving indexes. The purpose is to effectively clarify linkage relationship between indexes. A regression analysis method is used to quantify the mutual influence relationship. A multi-index coupling peak shaving optimization model for generation scheduling is developed and finally solved by using genetic algorithm. The invention can objectively describe comprehensive peak shaving requirements of a provincial power grid, and is helpful for generation dispatcher to reduce subjective assumptions in the decision-making. Compared with a conventional single-index peak-shaving model, peak-shaving effect of our method is obviously improved, which provides a novel technical way for peak-shaving generation scheduling and decision-making in short-term operations of power systems.

TECHNICAL FIELD

The invention relates to the field of hydropower scheduling, focusing ona multi-objective optimization method for power system peak-shavingbased on index linkage analysis.

BACKGROUND

Peak-shaving is a core task of the short-time generation scheduling ofpower systems, and it's also one of the key issues at home and abroad.Due to the large load difference between peak and off-peak and theincreasing intermittent energy grid connection, most provincial powergrids in China are facing severe peak-shaving pressure. Peak-shaving hasbecome one of the main targets in the daily power generation schedulingof these power grids. From the perspective of the power grid, the mainpurpose of peak-shaving is optimizing the power generation plans of thedispatchable power plants to reduce the peak-valley difference ofresidual loads and smooth the fluctuations in the net load curve forcoal-fired units with poor regulating performance. In order to achievethis goal, how to accurately describe the peak-shaving requirementsbecomes a key part of power generation scheduling.

Nowadays, the technical measures to deal with the peak-shaving pressureof power systems mainly include two aspects: One is to focus on theconstruction of energy storage systems and its potential ofpeak-shaving, such as introducing a peak demand charge into the energycost function, the joint operation optimization of a hybrid systemincluding a photovoltaic farm and battery-based energy storage, thecomplementary operation of photovoltaic-wind plant and hydro-pumpingstorage. The other is the peak-shaving optimization of conventionalpower plants such as thermal power plants, hydropower plants and so on.This kind of methods generally formulates models based on the variableof the net load. Some forms of a function of the net load are used todescribe the peak-shaving needs, such as the minimum of the maximumresidual peak-valley difference, the minimum of the variance or meansquare deviation of remaining load, the maximum of the peak-shavingcapacity and so on. Overall, these modelling methods can effectivelyreduce the maximum peak-valley difference or smooth the fluctuations inthe load curve in appropriate situations. However, for the peak-shavingpressure or requirements facing provincial power grids, it is usuallydescribed by a peak-shaving index such as the peak-valley difference,the mean square deviation of remaining load, etc. The comprehensiveanalysis of the load demand is lacked. When the load demand of differentpower girds changed at different dates, there is a lack of attention onthe applicability of each peak-shaving index and which peak-shavingindex should be adopted to optimize the generation schedules for betterpeak-shaving results. It's a very important problem in the research ofpeak operations.

For the above problem, the invention proposes a multi-objectiveoptimization method for power system peak-shaving based on index linkageanalysis. This work is supported by the National Natural ScienceFoundation of China (52079014). Based on the engineering background ofXiluodu hydropower plant transmitting power to ZheJiang Power Grid, theinvention is tested and applied. The results show that the invention caneffectively avoid the deviation of peak-shaving results caused bysingle-index optimization as well as the uncertainty and irrationalityof subjective factors. Thus, the optimized generation schedules are moreeffective and practical.

SUMMARY

The technical problem to be solved by the invention is to provide amulti-objective optimization method for power system peak-shaving basedon index linkage analysis. The purpose is to reduce uncertainty andirrationality caused by the weight coefficient selection in themulti-objective processing strategy and to obtain balanced peak-shavingoptimization results.

The technical solution of the invention:

A multi-objective optimization method for power system peak-shavingbased on index linkage analysis is characterized by the following steps:

(1) Set initial calculation conditions, including operation conditionsand constraints of hydropower plants, as well as conditions for electricand hydraulic operation demands;

(2) Five peak-shaving indexes are introduced, including load variance

${F_{1} = {\sum\limits_{i = 1}^{n}\left( {P_{r,i} - {\overset{¯}{P}}_{r}} \right)^{2}}},$load difference between peak and off-peak F₂=P_(r,max)−P_(r,min),maximum variation rate of loads F₃=|P_(r,i)−P_(r,i-1)|_(max), maximumload F₄=P_(r,max), minimum load F₅=P_(r,min); where P_(r) is residualload series; P_(r,i) is residual load in period i; P_(r,max) is maximumresidual load; P_(r,min) is minimum residual load; n is total number ofperiods. Thus, a multi-objective function is obtained:

$\begin{matrix}{{MOP}\left\{ \begin{matrix}{\min F_{1}} \\{\min\; F_{2}} \\{\min\; F_{3}} \\{\min\; F_{4}} \\{\max F_{5}}\end{matrix} \right.} & (1)\end{matrix}$

(3) According to historical data of daily load curve, values of theabove five indexes are calculated to obtain five time series that arenamed l₁, l₂, l₃, l₄ and l₅, respectively.

(4) The following formula is used to calculate correlation coefficientsbetween pairs of five series mentioned above:

$\begin{matrix}{r = \frac{\sum\limits_{i = 1}^{n}{\left( {l_{x,i} - \overset{\_}{l_{x}}} \right)\left( {l_{y,i} - \overset{\_}{l_{y}}} \right)}}{\sqrt{\sum\limits_{i = 1}^{n}\left( {l_{x,i} - \overset{\_}{l_{x}}} \right)^{2}}\sqrt{\sum\limits_{i = 1}^{n}\left( {l_{y,i} - \overset{¯}{l_{y}}} \right)^{2}}}} & (2)\end{matrix}$where r is correlation coefficient; l_(x,i) and l_(y,i) are values ofsequence l_(x) and sequence l_(y) in period i, respectively; l_(x) andl_(y) are average values of sequence l_(x) and sequence l_(y),respectively. Correlation between five indexes is judged according tothe follows: two indexes are highly correlated if 0.8<|r|<1; two indexesare moderately correlated if 0.3<|r|<0.8 and two indexes are lowlycorrelated if |r|<0.3.

(5) Based on the correlation between indexes, correlation analysis iscarried out for indexes with moderate or higher correlation; a weightcoefficient method is adopted for indexes with low correlation; amulti-objective peak-shaving objective is transformed into asingle-objective peak-shaving model; specific steps are as follows:

Step 1. Check stationarity of time series. ADF test method is used; ageneration process of sequence l is one of the following three forms:

$\begin{matrix}{{{\Delta\; l_{t}} = {{\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (3) \\{{{\Delta\; l_{t}} = {\mu + {\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (4) \\{{{\Delta l_{t}} = {\mu + {\alpha t} + {\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (5)\end{matrix}$where αt is time trend; μ is displacement item;

$\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta\; l_{t - 1}}$is lag item; k is number of lag items; Δl_(t) is increment of t insequence l that is calculated by Δl_(t)=l_(t)−l_(t−1); originalhypothesis and alternative hypothesis are

$\left\{ {\begin{matrix}{{H_{0}:\beta} = {0\left( {y_{t}\mspace{14mu}{nonstationary}} \right)}} \\{H_{1}:{\beta < {0\left( {y_{t}\mspace{14mu}{stationary}} \right)}}}\end{matrix},} \right.$respectively; a judging criteria is

$\left\{ {\begin{matrix}{{{{ADF}\mspace{14mu}{value}} \geq {threshold}},\mspace{14mu}{{accept}\mspace{14mu} H_{0}},{y_{t}\mspace{20mu}{nonstationary}}} \\{{{{ADF}\mspace{14mu}{value}} < {threshold}},\mspace{14mu}{{accept}\mspace{14mu} H_{0}},{y_{t}\mspace{14mu}{stationary}}}\end{matrix};} \right.$test the above three forms in turn; if one of them rejects originalhypothesis, the sequence would be stable;

Step 2. Granger causality test is conducted among sequences withmoderate correlation and stability. Regression models between them areas follows:

$\begin{matrix}{l_{y,t} = {{\sum\limits_{\;^{i = 1}}^{p}{\alpha_{i}l_{x,{t - i}}}} + {\sum\limits_{\;^{i = 1}}^{p}{\beta l_{y,{t - i}}}} + ɛ_{1t}}} & (6) \\{l_{x,t} = {{\sum\limits_{i = 1}^{p}{\lambda_{i}l_{y,{t - i}}}} + {\sum\limits_{i = 1}^{p}{\delta_{i}l_{x,{t - i}}}} + ɛ_{2t}}} & (7)\end{matrix}$where α_(i), β_(i), γ_(i) and δ_(i) are regression coefficients; ε_(1t),ε_(2t) are random error items; p is a displacement item; carry outhypothesis testing for the two models, respectively; original hypothesisand alternative hypothesis are

$\left\{ \begin{matrix}{{H_{0}:\alpha_{i}} = {0\mspace{11mu}\left( {{there}\mspace{14mu}{is}\mspace{14mu}{no}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}} \\{H_{1}:{\alpha_{i} \neq {0\mspace{11mu}\left( {{there}\mspace{14mu}{is}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}}}\end{matrix}\quad \right.\quad$and

$\left\{ \begin{matrix}{{H_{0}:\lambda_{i}} = {0\mspace{11mu}\left( {{there}\mspace{14mu}{is}\mspace{14mu}{no}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}} \\{H_{1}:{\lambda_{i} \neq {0\mspace{11mu}\left( {{there}\mspace{14mu}{is}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}}}\end{matrix}\quad \right.$respectively; according to test results, the Granger causality betweensequences could be obtained; if there is Granger causality between l_(x)and l_(y), then corresponding index X would has an impact on index Y.

Step 3. Co-integration of unstable sequences. If an unstable sequencel_(y) becomes a stable sequence after D times difference, then l_(y) isintegrated of order D; use OLS co-integration regression equation fortwo unstable sequences l_(x) and l_(y) with same integrated order; whereT, U are regression coefficients; X_(t) is residual; when X_(t) isstable, then corresponding index X would has an impact on index Y;

Step 4. Regression analysis based on Granger causality of stablesequence and quantitative variation relationship of unstable sequences;choose the most affected one of five indexes as θ₀; indexes which affectθ₀ are denoted as θ₁, θ₂, . . . , θ_(m), respectively; where m is numberof indexes that affect θ₀; a formula θ₀=ƒ(θ₁, θ₂, . . . , θ_(m)) iscalculated by regression analysis;

Step 5. Construct a peak-shaving objective function; based on formulaθ₀=ƒ(θ₁, θ₂, . . . , θ_(m)), the multi-objective function mentionedabove can be transformed into a single-objective function by using acoefficient weight method for residual indexes such as θ_(n+1), θ_(n+2),. . . , θ₄; the function expression is as follows:

$\begin{matrix}{{\min F} = {{f\left( {\theta_{1},\theta_{2},\ldots\mspace{14mu},\theta_{m}} \right)} + {\sum\limits_{\;^{i = {m + 1}}}^{4}{\beta_{i}\theta_{i}}}}} & (8)\end{matrix}$

(5) Get the generation profile of power plants by using GeneticAlgorithm method to solve the objective function in (4).

The invention has the following beneficial effects: compared with anormal single-objective peak-shaving model, a multi-objective couplingpeak-shaving models could quantitatively describe reasonablerelationship among five indexes of fluctuation degree, fluctuationrange, peak value, valley value and maximum variation rate of load,balance peak-shaving results; it could effectively avoid an error ofpeak-shaving caused by a single-objective optimization; differentindexes essentially reflect different peak-shaving demands of a powergrid; through a correlation analysis between indexes, it could specifyand definite load regulation demands of different power grids anddifferent times, and obtain a more reasonable peak-shaving model, makingoptimized generation schedules more effective and useful; the inventiondescribes complicated peak-shaving demand of power systems by usingeconomic theory methods; this research approach is benefit to promotedescription of practical problems in power system generation scheduling,which could provide relatively objective optimization results andeffectively avoid uncertainty and irrationality caused by subjectivefactors.

DESCRIPTION OF DRAWINGS

FIG. 1: A general solution frame diagram of the invention.

FIG. 2: A causal diagram of indexes.

FIG. 3: Analysis of fitting effect.

DETAILED DESCRIPTION

The specific embodiments of the present invention are further describedbelow in conjunction with the drawings and technical solutions.

In order to describe peak-shaving effect accurately, the inventioncomprehensively considers variance, peak-valley difference, maximum,minimum and maximum variation rate of residual loads. Among them, thevariance mainly reflects frequency of load changes in rise and falldirections, which is a requirement for flexible adjustment ability ofgenerating units:

$\begin{matrix}{\sigma^{2} = {\sum\limits_{i = 1}^{n}\left( {P_{r,i} - {\overset{\_}{P}}_{r}} \right)^{2}}} & (9)\end{matrix}$P _(r,i) =P ₁ −N ₁  (10)

where P_(r,i) is residual load at period i; P_(i) is load ofreceiving-end power grid; N_(i) is generation of optimized power plant.

The peak-valley difference Δ is a variation range of load, the smallerthe range is, the smaller the peak-shaving capacity demand of generatingunits is:Δ=P _(r,max) −P _(r,min)  (11)where P_(r,max) is maximum residual load; P_(r,min) is minimum residualload.

The maximum load N_(r,max) reflects peak-shaving pressure of peak loadperiods that is a period with the maximum load response pressure ofgenerating units. To ensure safe and stable operations of power grids,the residual load is expected to stay at a low level during peak periodsin order to avoid peak-shaving difficulty of generating units.N _(max) =P _(r,max)  (12)

The minimum load N_(r,min) reflects peak-shaving pressure of valley loadperiods, which is limited by a forced generation of non-regulatingunits. During low load periods, the output of generating units isdifficult to be effectively reduced with load, and it forces thermalunits to work with in-depth peak-shaving, thereby resulting in asignificant increase in energy consumption. Therefore, the residual loadis expected to stay at a high level during valley periods.N _(min) =P _(r,min)  (13)

The maximum variation rate ΔP_(r,max) of load is maximum variation ofload in adjacent periods that reflects demand of load for climbingability of generating units. The smaller the maximum variation rate ofload is, the lower the demand for climbing ability of generating unitsis.ΔP _(r,max) =|P _(r,i) −P _(r,i-1)|_(max)  (14)

The objective function is constructed to minimize the variance,peak-valley difference, maximum and maximum variation rate of residualloads, and maximize the minimum of residual loads to ensure peak-shavingdemands of different loads.

$\begin{matrix}{{MOP}\left\{ \begin{matrix}{\min\sigma^{2}} \\{\min\Delta} \\{\min N_{\max}} \\{\max N_{\min}} \\{\min\Delta P_{r,\max}}\end{matrix} \right.} & (15)\end{matrix}$

The objective function mentioned above should meet the followingconstraints:

Water balance equation:V _(t+1) =V _(t)+3600(Q _(t) −S _(t))Δt  (16)

-   -   where S_(t)=q_(t)+d_(t);    -   V_(t) is reservoir capacity of power plant at time t; Q_(t) is        reservoir inflow in period t, m³/s; S_(t) is storage outflow of        power plant in period t, m³/s; q_(t) is power discharge of power        plant in period t, m³/s; d_(t) is abandoning water of power        plant in period t, m³/s.

Relationship of water level and reservoir capacity:Z _(t)=ƒ_(ZV)(V _(t))  (17)where Z_(t) is water level of power plant in period t, m; ƒ_(ZV)(⋅) is arelation curve of water level and reservoir capacity of power plant.

The relationship of tailrace level and discharge capacity:Z _(t) ^(d)=ƒ_(z) _(d) (S _(t))  (18)where Z_(t) ^(d) is tailrace level of power plant in period t, m; ƒ_(Z)_(d) (⋅) is a relation curve of tailrace level and discharge capacity ofpower plant.

Water head relationship:h _(t)=(Z _(t) +Z _(t+1))/2−Z _(t) ^(d)  (19)where h_(t) is water head of power plant in period t, m.

The relationship of power station output:N _(t)=ƒ_(qh)(q _(t) ,h _(t))  (20)where N_(t) is generation of power plant in period t, MW; ƒ_(qh)(⋅) is aNQH curve of power plant.

Electricity generation control:

$\begin{matrix}{E_{total} = {\sum\limits_{t = 1}^{T}{N_{t} \times \Delta t}}} & (21)\end{matrix}$where E_(total) is given total amount of power generation, MWh.

Constraint of reservoir water level:Z _(t) ≤Z _(t) ≤Z _(t)  (22)where Z _(t) and Z _(t) are upper and lower water level of reservoir inperiod t, respectively, m.

Upper and lower limits of power generation:N _(t) ≤N _(t) ≤N _(t)  (23)where N _(t) and N _(t) are upper and lower average generation in periodt of power plant, respectively, MW.

Constraint of power discharge:q _(t) ≤q _(t) ≤q _(t)  (24)where q _(t) and q _(t) are upper and lower power discharge in period tof power plant, respectively, m³/s.

Constraint of storage outflow:S _(t) ≤S _(t) ≤S _(t)  (25)where S _(t) and S _(t) are upper and lower storage outflow in period tof m power plant, m³/s.

Taking the above multi-objective optimization model as the background,the invention is applied in a complete way according to the following(1) to (5) steps:

(1) Set initial calculation conditions, including operation conditionsand constraints of hydropower plants, as well as conditions for electricand hydraulic operation demands;

(2) Five peak-shaving indexes are introduced, including load variance

${F_{1} = {\sum\limits_{i = 1}^{n}\left( {P_{r,i} - {\overset{\_}{P}}_{r}} \right)^{2}}},$load difference between peak and off-peak F₂=P_(r,max)−P_(r,min),maximum variation rate of loads F₃=|P_(r,i)−P_(r,i-1)|_(max), maximumload F₄=P_(r,max), minimum load F₅=P_(r,min); where P_(r) is residualload series; P_(r,i) is residual load in period i; P_(r,max) is maximumresidual load; P_(r,min), is minimum residual load; n is total number ofperiods. Thus, we get a multi-objective function:

$\begin{matrix}{{MOP}\left\{ \begin{matrix}{\min\; F_{1}} \\{\min\; F_{2}} \\{\min F_{3}} \\{\min\; F_{4}} \\{\max F_{5}}\end{matrix} \right.} & (26)\end{matrix}$

(3) According to historical data of daily load curve, we can calculatevalues of the above five indexes to obtain five time series that arenamed l₁, l₂, l₃, l₄ and l₅, respectively.

(4) The following formula is used to calculate correlation coefficientsbetween pairs of five series mentioned above:

$\begin{matrix}{r = \frac{\sum\limits_{i = 1}^{n}{\left( {l_{x,i} - \overset{\_}{l_{x}}} \right)\left( {l_{y,i} - \overset{\_}{l_{y}}} \right)}}{\sqrt{\sum\limits_{i = 1}^{n}\left( {l_{x,i} - \overset{\_}{l_{x}}} \right)^{2}}\sqrt{\sum\limits_{i = 1}^{n}\left( {l_{y,i} - \overset{\_}{l_{y}}} \right)^{2}}}} & (27)\end{matrix}$where r is correlation coefficient; l_(x,i) and l_(y,i) are values ofsequence l_(x) and sequence l_(y) in period i, respectively; l_(x) andl_(y) are average values of sequence l_(x) and sequence l_(y),respectively. Correlation between five indexes is judged according tothe follows: two indexes are highly correlated if 0.8<|r|<1; two indexesare moderately correlated if 0.3<|r|<0.8 and two indexes are lowlycorrelated if |r|<0.3.

(5) Based on the correlation between indexes, correlation analysis iscarried out for indexes with moderate or higher correlation; a weightcoefficient method is adopted for indexes with low correlation; amulti-objective peak-shaving objective is transformed into asingle-objective peak-shaving model; specific steps are as follows:

Step 1. Check stationarity of time series. ADF test method is used; ageneration process of sequence l is one of the following three forms:

$\begin{matrix}{{{\Delta\; l_{t}} = {{\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta\; l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (28) \\{{{\Delta\; l_{t}} = {\mu + {\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta\; l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (29) \\{{{\Delta\; l_{t}} = {\mu + {\alpha t} + {\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta\; l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (30)\end{matrix}$where αt is time trend; μ is displacement item;

$\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta\; l_{t - i}}$is lag item; k is number of lag items; Δl_(t) is increment of t insequence l that is calculated by Δl_(t)=l_(t)−l_(t−1); originalhypothesis and alternative hypothesis are

$\left\{ {\begin{matrix}{{H_{0}:\beta} = {0\left( {y_{t}\mspace{14mu}{nonstationary}} \right)}} \\{H_{1}:{\beta < {0\left( {y_{t}\mspace{14mu}{stationary}} \right)}}}\end{matrix},}\quad \right.$respectively; a judging criteria is

$\left\{ {\begin{matrix}{{{{ADF}\mspace{14mu}{value}} \geq {threshold}},\mspace{14mu}{{accept}\mspace{14mu} H_{0}},{y_{t}\mspace{20mu}{nonstationary}}} \\{{{{ADF}\mspace{14mu}{value}} < {threshold}},\mspace{14mu}{{accept}\mspace{14mu} H_{0}},{y_{t}\mspace{14mu}{stationary}}}\end{matrix};}\quad \right.$test the above three forms in turn; if one of them rejects originalhypothesis, the sequence would be stable;

Step 2. Granger causality test is conducted among sequences withmoderate correlation and stability. Regression models between them areas follows:

$\begin{matrix}{l_{y,t} = {{\sum\limits_{i = 1}^{p}{\alpha_{i}l_{x,{t - i}}}} + {\sum\limits_{i = 1}^{p}{\beta_{i}l_{y,{t - i}}}} + ɛ_{1t}}} & (31) \\{l_{x,t} = {{\sum\limits_{i = 1}^{p}{\lambda_{i}l_{y,{t - 1}}}} + {\sum\limits_{i = 1}^{p}{\delta_{i}l_{x,{t - i}}}} + ɛ_{2t}}} & (32)\end{matrix}$where α_(i), β_(i), γ_(i) and δ_(i) are regression coefficients; ε_(1t),ε_(2t) are random error items; p is a displacement item; carry outhypothesis testing for the two models, respectively; original hypothesisand alternative hypothesis are

$\left\{ \begin{matrix}{{H_{0}:\alpha_{i}} = {0\mspace{11mu}\left( {{there}\mspace{14mu}{is}\mspace{14mu}{no}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}} \\{H_{1}:{\alpha_{i} \neq {0\mspace{11mu}\left( {{there}\mspace{14mu}{is}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}}}\end{matrix}\quad \right.$and

$\left\{ \begin{matrix}{{H_{0}:\lambda_{i}} = {0\mspace{11mu}\left( {{there}\mspace{14mu}{is}\mspace{14mu}{no}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}} \\{H_{1}:{\lambda_{i} \neq {0\mspace{11mu}\left( {{there}\mspace{14mu}{is}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}}}\end{matrix}\quad \right.$respectively; according to test results, the Granger causality betweensequences could be obtained; if there is Granger causality between l_(x)and l_(y), then corresponding index X would has an impact on index Y.

Step 3. Co-integration of unstable sequences. If an unstable sequencel_(y) becomes a stable sequence after D times difference, then l_(y) isintegrated of order D; use OLS co-integration regression equation fortwo unstable sequences l_(x) and l_(y) with same integrated order; whereT, U are regression coefficients; X_(t) is residual; when X_(t) isstable, then corresponding index X would has an impact on index Y;

Step 4. Regression analysis based on Granger causality of stablesequence and quantitative variation relationship of unstable sequences;choose the most affected one of five indexes as θ₀; indexes which affectθ₀ are denoted as θ₁, θ₂, . . . , θ_(m), respectively; where m is numberof indexes that affect θ₀; a formula θ₀=ƒ(θ₁, θ₂, . . . , θ_(m)) iscalculated by regression analysis;

Step 5. Construct a peak-shaving objective function; based on formulaθ₀=ƒ(θ₁, θ₂, . . . , θ_(m)), the multi-objective function mentionedabove can be transformed into a single-objective function by using acoefficient weight method for residual indexes such as θ_(n+1), θ_(n+2),. . . , θ₄; the function expression is as follows:

$\begin{matrix}{{\min F} = {{f\left( {\theta_{1},\theta_{2},\ldots\mspace{14mu},\theta_{m}} \right)} + {\sum\limits_{i = {m + 1}}^{4}{\beta_{i}\theta_{i}}}}} & (33)\end{matrix}$

(5) Get the generation profile of power plants by using GeneticAlgorithm method to solve the objective function in (4).

Now, the presented model was verified by applying it to an actualengineering of Xiluodu Hydropower Plant. Xiluodu is one of main powersuppliers of “Electricity Transmission Project from West to East”, andit's also the biggest hydropower plant on the Jinsha River. There are 9identical hydropower units of 770 MW each installed on the left andright sides of the river. Among them, generating units of the plant onthe left banks provide electricity to Zhejiang via the Binjin UHVDCline. The proportion of thermal power in the Zhejiang Power Grid is upto 95%, which lacks flexible generation capacity and exhibits largedaily load differences between peak and off-peak periods, facing a highpeak-shaving pressure. Therefore, it's necessary to make use of flexibleregulation characteristics of hydropower to meet peak-shavingrequirements. Typical daily load curves of summer and winter in Zhejiangprovince are selected for calculation. The calculation results ofcorrelation coefficient between indexes are given in Table 2, amongwhich the maximum variation rate is of low correlation with otherindexes, thus only four indexes, variance, peak-valley difference,maximum, minimum, need to be analyzed jointly. Stability test resultsare listed in Table 3, and it can be obtained that four indexes arestable without co-integration. The results of Granger causality test arelisted in Table 4, so a causal guidance diagram can be drawn betweenthese indexes, as shown in FIG. 2. Regression analysis was conductedaccording to the causality, and the following fitting formula could beobtained, and the fitting effect is shown as FIG. 3:θ₀=9.23×10⁸−0.041×θ₁ ²+0.462×θ₂ ²−37535×θ₂+0.132×θ₃ ²−9529×θ₃where θ₀ is variance; θ₁ is peak-valley difference; θ₂ is minimum; θ₃ ismaximum.

The index of maximum variation rate is treated by a weight coefficientmethod, and the following function could be obtained:min F=9.23×10⁸−0.041×θ₁ ²+0.462×θ₂ ²−37535×θ₂+0.132×θ₃ ²−9529×θ₃+βθ₄where θ₄ is maximum variation rate; β is weight coefficient.

Using the objective function mentioned above for the optimization, theindexes and decreasing ranges of corresponding index could be obtainedand listed in Table 5 and Table 6, respectively. Among them, indexes ofvariance, peak-valley difference, maximum and the maximum variation rateshowed a significant decline and the index of minimum showed a smalldecline. It can be seen that it has a significant peak-shaving effect insummer and winter, and effectively alleviates the peak-shaving pressureof the power grid.

TABLE 1 Basic data of power plant Normal Maximum Maximum Plant Deadwater high water Installed Regulation reservoir power name level/mlevel/m capacity/MW performance discharge/m³/s discharge/m³/s Xiluodu374.00 600.00 13860 multi-year 100000 7749 regulating

TABLE 2 Calculation results of correlation coefficient MaximumFluctuation Fluctuation Peak Valley variation R degree range value valuerate Variance 1 Peak-valley 0.949 1 difference Maximum 0.720 0.861 1Minimum −0.950 −0.922 −0.597 1 Maximum 0.140 0.276 0.094 −0.056 1variation rate

TABLE 3 Stability testing ADF test 1% critical 5% critical value valuevalue C1onclusion Variance −4.42021 −3.43927 −2.8653 stable Peak-valleydifference −3.74418 −3.43926 −2.86536 stable Maximum −4.4423 −3.43926−2.86536 stable Minimum −4.13918 −3.43926 −2.86536 stable

TABLE 4 Results of Granger F- Null hypothesis statistic P valuePeak-valley difference is not Granger cause 25.349 5  3.00 × 10−27 ofvariance Variance is not Granger cause of Peak- 13.853 4  6.00 × 10−15valley difference Minimum is not Granger cause of maximum 2.335 3 0.0306 Variance is not Granger cause of minimum 1.616 2 0.139 8 Maximum isnot Granger cause of variance 18.092 4  1.00 × 10−19 Variance is notGranger cause of maximum 2.499 2 0.021 2 Maximum is not Granger cause ofPeak- 4.781 8 9.00 × 10−05 valley difference Peak-valley difference isnot Granger cause 0.488 7 0.817 1 of maximum Minimum is not Grangercause of Peak- 4.781 8 9.00 × 10−05 valley difference Peak-valleydifference is not Granger cause 14.320 0  2.00 × 10−15 of minimumMinimum is not Granger cause of maximum 0.488 7 0.817 1 Maximum is notGranger cause of minimum 14.320 0  2.00 × 10−15

TABLE 5 Calculation results of indexes Index θ₀/MW² θ₁/MW θ₂/MW θ₃/MWθ₄/MW Typical Original 17 612 306  13 262  37 982 51 244 4 696 summerload days Residual 4 520 608 6 333 37 982 44 315 1 616 load TypicalOriginal 2 531 061 5 156 43 401 48 557 4 183 winter load days Residual  144 675 1 103 41 789 42 892 1 103 load

TABLE 6 Decreasing range of indexes (%) Index θ₀ θ₁ θ₂ θ₃ θ₄ Typicalsummer days 74.3 52.2 0.0 13.5 65.6 Typical winter days 94.3 78.6 3.711.7 73.6

The invention claimed is:
 1. A multi-objective optimization method forpower system peak-shaving based on index linkage analysis, whereincomprising the following steps: (1) set initial calculation conditions,including operation conditions and constraints of hydropower plants, aswell as conditions for electric and hydraulic operation demands; (2)five peak-shaving indexes are introduced, including load variance${F_{1} = {\sum\limits_{1 = 1}^{n}\left( {P_{r,i} - {\overset{\_}{P}}_{r}} \right)^{2}}},$load difference between peak and off-peak F₂=P_(r,max)−P_(r,min),maximum variation rate of loads F₃=|P_(r,i)−P_(r,i-1)|_(max), maximumload F₄=P_(r,max), minimum load F₅=P_(r,min); where P_(r) is residualload series; P_(r,i) is residual load in period i; P_(r,max) is maximumresidual load; P_(r,min) is minimum residual load; n is total number ofperiods; thus, a multi-objective function is obtained: $\begin{matrix}{{MOP}\left\{ \begin{matrix}{\min\; F_{1}} \\{\min\; F_{2}} \\{\min F_{3}} \\{\min\; F_{4}} \\{\max F_{5}}\end{matrix} \right.} & (1)\end{matrix}$ (3) according to historical data of daily load curve,values of the above five indexes are calculated to obtain five timeseries that are named l₁, l₂, l₃, l₄ and l₅, respectively; (4) Eq. (2)is used to calculate correlation coefficients between pairs of fiveseries mentioned above: $\begin{matrix}{r = \frac{\sum\limits_{i = 1}^{n}{\left( {l_{x,i} - \overset{\_}{l_{x}}} \right)\left( {l_{y,i} - \overset{\_}{l_{y}}} \right)}}{\sqrt{\sum\limits_{i = 1}^{n}\left( {l_{x,i} - \overset{\_}{l_{x}}} \right)^{2}}\sqrt{\sum\limits_{i = 1}^{n}\left( {l_{y,i} - \overset{\_}{l_{y}}} \right)^{2}}}} & (2)\end{matrix}$ where r is correlation coefficient; l_(x,i) and l_(y,i)are values of sequence l_(x) and sequence l_(y) in period i,respectively; l_(x) and l_(y) are average values of sequence l_(x) andsequence l_(y), respectively; correlation between five indexes is judgedaccording to the follows: two indexes are highly correlated if0.8<|r|<1; two indexes are moderately correlated if 0.3<|r|<0.8 and twoindexes are lowly correlated if |r|<0.3; (5) based on the correlationbetween indexes, correlation analysis is carried out for indexes withmoderate or higher correlation; a weight coefficient method is adoptedfor indexes with low correlation; a multi-objective peak-shavingobjective is transformed into a single-objective peak-shaving model;specific steps are as follows: Step
 1. check stationarity of timeseries; ADF test method is used; a generation process of sequence l isone of the following three forms: $\begin{matrix}{{{\Delta\; l_{t}} = {{\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta\; l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (3) \\{{{\Delta\; l_{t}} = {\mu + {\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta\; l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (4) \\{{{\Delta\; l_{t}} = {\mu + {\alpha t} + {\beta l_{t - 1}} + {\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta\; l_{t - i}}} + u_{i}}},{l_{0} = 0}} & (5)\end{matrix}$ where αt is time trend; μ is displacement item;$\sum\limits_{i = 1}^{k}{\gamma_{i}\Delta\; l_{t - i}}$ is lag item; kis number of lag items; Δl_(t) is increment of t in sequence l that iscalculated by Δl_(t)=l_(t)−l_(t−1); original hypothesis and alternativehypothesis are $\left\{ {\begin{matrix}{{H_{0}:\beta} = {0\left( {y_{t}\mspace{14mu}{nonstationary}} \right)}} \\{H_{1}:{\beta < {0\left( {y_{t}\mspace{14mu}{stationary}} \right)}}}\end{matrix},} \right.$ respectively; a judging criteria is$\left\{ {\begin{matrix}{{{{ADF}\mspace{14mu}{value}} \geq {threshold}},\mspace{14mu}{{accept}\mspace{14mu} H_{0}},{y_{t}\mspace{20mu}{nonstationary}}} \\{{{{ADF}\mspace{14mu}{value}} < {threshold}},\mspace{14mu}{{accept}\mspace{14mu} H_{0}},{y_{t}\mspace{14mu}{stationary}}}\end{matrix};} \right.$ test the above three forms in turn; if one ofthem rejects original hypothesis, the sequence would be stable; Step 2.Granger causality test is conducted among sequences with moderatecorrelation and stability; regression models between them are asfollows: $\begin{matrix}{l_{y,t} = {{\sum\limits_{i = 1}^{p}\;{\alpha_{i}l_{x,{t - i}}}} + {\sum\limits_{i = 1}^{p}\;{\beta_{i}l_{y,{t - i}}}} + ɛ_{1t}}} & (6)\end{matrix}$ $\begin{matrix}{l_{x,t} = {{\sum\limits_{i = 1}^{p}\;{\lambda_{i}l_{y,{t - i}}}} + {\sum\limits_{i = 1}^{p}\;{\delta_{i}l_{x,{t - i}}}} + ɛ_{2t}}} & (7)\end{matrix}$ where α_(i), β_(i), γ_(i) and δ_(i) are regressioncoefficients; ε_(1t), ε_(2t) are random error items; p is a displacementitem; carry out hypothesis testing for the two models, respectively;original hypothesis and alternative hypothesis are$\left\{ \begin{matrix}{{H_{0}:\alpha_{i}} = {0\mspace{11mu}\left( {{there}\mspace{14mu}{is}\mspace{14mu}{no}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}} \\{H_{1}:{\alpha_{i} \neq {0\mspace{11mu}\left( {{there}\mspace{14mu}{is}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}}}\end{matrix}\quad \right.$ and $\left\{ \begin{matrix}{{H_{0}:\lambda_{i}} = {0\mspace{11mu}\left( {{there}\mspace{14mu}{is}\mspace{14mu}{no}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}} \\{H_{1}:{\lambda_{i} \neq {0\mspace{11mu}\left( {{there}\mspace{14mu}{is}\mspace{14mu}{Granger}\mspace{14mu}{causality}\mspace{14mu}{relation}\mspace{14mu}{for}\mspace{14mu} l_{x}\mspace{14mu}{by}\mspace{14mu} l_{y}} \right)}}}\end{matrix}\quad \right.$ respectively; according to test results, theGranger causality between sequences could be obtained; if there isGranger causality between l_(x) and l_(y), then corresponding index Xwould has an impact on index Y; Step
 3. co-integration of unstablesequences; if an unstable sequence l_(y) becomes a stable sequence afterD times difference, then l_(y) is integrated of order D; use OLSco-integration regression equation for two unstable sequences l_(x) andl_(y) with same integrated order; where T, U are regressioncoefficients; X_(t) is residual; when X_(t) is stable, thencorresponding index X would has an impact on index Y; Step
 4. regressionanalysis based on Granger causality of stable sequence and quantitativevariation relationship of unstable sequences; choose the most affectedone of five indexes as θ₀; indexes which affect θ₀ are denoted as θ₁,θ₂, . . . , θ_(m), respectively; where m is number of indexes thataffect θ₀; a formula θ₀=ƒ(θ₁, θ₂, . . . , θ_(m)) is calculated byregression analysis; Step
 5. construct a peak-shaving objectivefunction; based on Eq. θ₀=ƒ(θ₁, θ₂, . . . , θ_(m)), the multi-objectivefunction mentioned above can be transformed into a single-objectivefunction by using a coefficient weight method for residual indexes suchas θ_(n+1), θ_(n+2), . . . , θ₄; it is shown as follows: $\begin{matrix}{{\min\; F} = {{f\left( {\theta_{1},\theta_{2},\ldots,\theta_{m}} \right)} + {\sum\limits_{i = {m + 1}}^{4}\;{\beta_{i}\theta_{i}}}}} & (8)\end{matrix}$ (5) get the generation profile of power plants by usinggenetic algorithm method to solve the objective function in (4).